TD Dynamique adaptative - MDet

Modele ecologique

On considere un modèle de competition. Dans notre modèle, 21 phénotypes dont le trait X est distribué sur un intervalle allant de 0 à 2, sont en compétition. Chaque population de chaque phénotype possède un taux de croissance intrasèque identique r. La densité maximale est spécifique à chaque phénotype et dépend uniquement de la valeur de son trait selon la formule :

\[ K(x) = max(K_0 - λ(x-x_0)^2, 0) \] L’intensité de la compétition entre deux phénotypes i et j vaut : \[ a(x_i, x_j) =exp(-\frac{(x_i-xj)^2}{2σ^2}) \]

La dynamique de chaque population s’écrit alors :

\[ \frac{dn_i}{dt} = r.n_i.(\sum_{j=1}^{n} \frac{a(x_i, x_j).n_j}{K(x_i)} ) \] Pour l’ensemble des manipulations, on fixe certaines valeurs par défaut aux paramètres.

#Parametres
i = 21 # Nombre d'especes
Xmin = 0
Xmax = 2
X0 = (Xmax - Xmin)/2
K0 = 1
lambda = 1
sigma = 1 #etalement de la competition
r = 1

On implèmente les fonctions mathématiques du modèle.

############ Fonctions #############

K_cap = function(x, K0_ = K0, lambda_= lambda, x0 = X0){
   (K0_ - lambda_*(x-x0)**2)
}

a = function(x1,x2, sigma_a = sigma){ 
  exp( (-1/2*(x2-x1)**2)/sigma_a**2 )
}

fitness = function(x1, x2, r_ = r ){
  r_*(1-a(x1, x2)*(K_cap(x1)/K_cap(x2)))
}

Fitness d’invasion

Afin de comprendre ce qu’implique notre modèle avec les valeurs de ses paramètres initiaux, nous représentons graphiquement la valeur de la fitness relative dans le cas de l’apparition d’un nouveau phénotype. Ces figure en trois dimensions fournissent le paysage adaptatif de notre modèle.

## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 25);
x2_interval  <- seq(0.1, 1.9, length.out = 25);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

p  <- persp(x1_interval,x2_interval,z, theta = 30, phi = 20,
            col = "lightblue", shade = 0.4, ticktype = "detailed")

plot_ly() %>% add_surface(x = ~x1_interval, y = ~x2_interval, z = ~z)

Le Pairwise Invasibily Plot (PIP) fournit un résumé de la situation en représentant seulement le signe de la différence des fitness.

## Set range and domain of plot
x1_interval  <- seq(0.1, 1.9, length.out = 500);
x2_interval  <- seq(0.1, 1.9, length.out = 500);

## Interpolate surface

z  <- outer(x1_interval,x2_interval,
            FUN = fitness)

couleurs <- ifelse(z > 0, "lightblue", "black")

matrice_couleurs = couleurs
df <- data.frame(x = rep(1:nrow(matrice_couleurs), each = ncol(matrice_couleurs)),
                 y = rep(1:ncol(matrice_couleurs), times = nrow(matrice_couleurs)),
                 couleur = as.vector(matrice_couleurs))

ggplot(df, aes(x = x, y = y, fill = couleur)) +
  geom_tile() +
  scale_fill_identity() +
  labs(title = "PIP") +
  theme_minimal()+
  labs(fill = "Légende\nNoir = Négatif\nBleu = Positif")

Simulations

Pour simuler ces compétitions, le modèle est implémenté sous forme matricielle. La dynamique est simulée pour un temps de 1000 et pour deux valeurs de sigma.

############### Simuler LV pour un nombre arbitraire d'espece ############

time_max = 1000

#Conditions initiales
x = seq(from = 0, to = 2, length.out = i) # Valeurs des phenotypes
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}# K de chaque espece

n = rep(1/i,i) #densite de pop initiale

# On met sous forme de matrice
N0 = t(as.matrix(n))
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)
X = matrix(rep(x,i), ncol = i, nrow = i, byrow = TRUE)

#Calcul de la matrice M
D = t(X) - X


#Resolution du systeme d'ODE
model <- function(t,N,sigma){
  
  M = exp(-0.5*D^2/(sigma^2))/t(K)
  
  dN <- r * N * (1 - N %*% M)
  
  return(list(dN))
}

ode <-ode(N0,0:time_max,model,parms = sigma)

sigma = 5
# Representation

ode <- ode(N0,0:time_max,model,parms = sigma)

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))

sigma = 0.3
# Representation

ode <-ode(N0,0:time_max,model,parms = sigma)

ode <- as.data.frame(ode)
#colnames(ode) <- c("t",'x1','x2','x3')
ode_plot = ode%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ggplot(ode_plot)+
  geom_line(aes(time, dens, col = sp_ID))

On observe que le nombre de phénotypes sui persistent dans le temps dépend notamment de la valeurs de sigma. Pour des faibles valeurs, un plus grand nombres de phénotypes se maintiennent.

ode_plot$sp_ID <- as.numeric(ode_plot$sp_ID)
ode_plot$trait <- x[as.numeric(ode_plot$sp_ID)]

ggplot(ode_plot)+
  geom_raster(aes(trait, time,  fill=dens**2))+
  scale_fill_gradient2(low = "white" ,
                       high = "red")

Le script et le graphique suivant permettent d’étudier le nombre de phénotypes se maintenant en fonction de la valeur de sigma. Les valeurs de sigma testées sont comprises entre 0.1 et 2.

seuil <- 0.001
sigma_bank = seq(0.1,2,by = 0.05)

nb_especes <- c()

for (sigma_i in sigma_bank){{
  
  ode2 <-ode(N0,0:time_max,model,parms = sigma_i)
  nb = sum(ode2[nrow((ode2)),-1] > seuil)
  }
  nb_especes <- c(c(nb_especes),c(nb))
  
}

res = tibble(sigma_bank,nb_especes)

ggplot()+
  geom_point(data = res, aes(x = sigma_bank, y = nb_especes))+
  labs(title = " ",
       x = "Sigma",
       y = "Nombre d'espèces en régime stationnaire")

Modelisation de mutations

Dans cette partie, on modélise la dynamique adaptative sur le temps long (T) d’une population en considérant que p = 10% de la population mute à intervalle de temps régulier (t). La population initiale est monomorphe et l’étalement de la compétition (sigma = 0.3) est assez faible.

#Implementation des mutations

#parametres
p = 0.1 #taux de la pop qui mute
Temps = 10000 #temps de la simulation
t = 400 #pas de temps entre chaque mutation
n = rep(0,i) 
n[4] <- 1 # densites de pop initiale (une seule espece)
N0 <- t(as.matrix(n))
sigma = 0.1

#Initialisation des valeures de K
k = c()
for (i in 1:length(x)){
  k = c(c(k),c(K_cap(x[i])))
}
K = matrix(rep(k,i), ncol = i, nrow = i, byrow = FALSE)

A = exp(-0.5*D^2/(sigma^2))
M = A/t(K)


# Evolution sans mutation sur le premier pas de temps t
ode <-ode(N0,1:t,model,parms = sigma, method = 'lsoda')
ode_mutation <- ode


seuil = 0.05 #densité min on considère qu'un trait avec une densité < seuil a une densité de 0


for (z in 2:(Temps/t)){
  
  #Elimination des traits sous le seuil
  tail <- as.matrix(tail(ode,1)[,1:i+1])
  for (j in 1:length(tail)){
    if(tail[j] < seuil){
      tail[j] <- 0
    }
  }
  N0 <- t(tail)
  
  # Mutation des traits restants
  for (index in which(tail != 0)){
    a = runif(1)
    if (a<0.5){ #mutation vers la gauche
      N0[max(0,index-1)] <- p*N0[index]
      N0[index] <- (1-p)*N0[index]
    } else { #mutation vers la droite
      N0[min(i,index+1)] <- p*N0[index]
      N0[index] <- (1-p)*N0[index]
    }
   
  }
  
  
  ode <-ode(N0,(z*t):(z*t+t),model,parms = sigma, method='lsoda')
  
  ode_mutation <- rbind(ode_mutation,ode) 

}

ode_mutation_plot = as.data.frame(ode_mutation)%>%
  pivot_longer(-time, values_to = "dens", names_to = "sp_ID")

ode_mutation_plot$sp_ID <- as.numeric(ode_mutation_plot$sp_ID)
ode_mutation_plot$trait <- x[as.numeric(ode_mutation_plot$sp_ID)]

 ggplot(ode_mutation_plot)+
   geom_tile(aes(trait, time, fill = dens))+
   scale_fill_gradient2(low = "white" ,
                        high = "red") +
   main_theme